#ifdef PYBINDINGS
#define NPY_NO_DEPRECATED_API  NPY_1_7_API_VERSION
#include <Python.h>
#include <numpy/arrayobject.h>
#endif

#include <stdio.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_odeiv2.h>

// outdated: http://www.scipy.org/Cookbook/C_Extensions/NumPy_arrays
// outdated: http://dsnra.jpl.nasa.gov/software/Python/numpydoc/numpy-13.html
// good: http://dan.iel.fm/posts/python-c-extensions/
// https://gist.github.com/dfm/3343607

int dyn (double t, const double y[], double f[], void *params);
int jac (double t, const double y[], double *dfdy, double dfdt[], void *params);

struct Param {
  gsl_odeiv2_driver *d;
  double t;
  double y[2];
};

struct Param p;



int dyn (double t, const double y[], double f[], void *params) {
  double mu = *(double *)params;
  f[0] = y[1];
  f[1] = -y[0] - mu*y[1]*(y[0]*y[0] - 1);
  return GSL_SUCCESS;
}

int jac (double t, const double y[], double *dfdy, double dfdt[], void *params) {
  double mu = *(double *)params;
  gsl_matrix_view dfdy_mat 
    = gsl_matrix_view_array (dfdy, 2, 2);
  gsl_matrix * m = &dfdy_mat.matrix; 
  gsl_matrix_set (m, 0, 0, 0.0);
  gsl_matrix_set (m, 0, 1, 1.0);
  gsl_matrix_set (m, 1, 0, -2.0*mu*y[0]*y[1] - 1.0);
  gsl_matrix_set (m, 1, 1, -mu*(y[0]*y[0] - 1.0));
  dfdt[0] = 0.0;
  dfdt[1] = 0.0;
  return GSL_SUCCESS;
}

#ifndef PYBINDINGS
int main(int argc, char** argv) {
  printf("hello\n");
  return 0;
}
#else

static PyObject* dynamic_model_init(PyObject* self, PyObject* args) {
  gsl_odeiv2_system sys = {dyn, jac, 2, &p};
  p.d =  gsl_odeiv2_driver_alloc_y_new(&sys, gsl_odeiv2_step_rk8pd, 1e-6, 1e-6, 0.0);
  p.t = 0.;
  p.y[0] = 1.; p.y[1] = 2.;
  return Py_BuildValue("i", 1);
}

static PyObject* dynamic_model_run(PyObject* self, PyObject* args) {
  PyObject *x_obj, *y_obj;
  double dt;
  if (!PyArg_ParseTuple(args, "OdO", &x_obj, &dt, &y_obj))
    return NULL;
  /* Interpret the input objects as numpy arrays. */
  PyObject *x_array = PyArray_FROM_OTF(x_obj, NPY_DOUBLE, NPY_ARRAY_IN_ARRAY);

#if 0
  /* How many data points are there? */
  int N = (int)PyArray_DIM(x_array, 0);

  /* Get pointers to the data as C-types. */
  double *x    = (double*)PyArray_DATA(x_array);
#endif
  //  for (int i=0; i<N; i++)
  //    printf("%f ", x[i]);
  //  printf("\n");

  //  x[0] =42.;

#if 0
  int status = gsl_odeiv2_driver_apply (p.d, &p.t, dt, p.y);
  if (status != GSL_SUCCESS){
    printf ("error, return value=%d\n", status);
  }
#endif
  /* Clean up. */
  Py_DECREF(x_array);
  return Py_BuildValue("i", 1);
}


/*
 * Bind Python function names to our C functions
 */
static PyMethodDef dynamic_model_methods[] = {
  {"init", dynamic_model_init, METH_VARARGS},
  {"run",  dynamic_model_run,  METH_VARARGS},
  {NULL, NULL}
};

/*
 * Python calls this to let us initialize our module
 */
void initdynamic_model_c() {
  (void) Py_InitModule("dynamic_model_c", dynamic_model_methods);
  /* Load `numpy` functionality. */
  import_array();
}

#endif /* PYBINDINGS */
